#########################################################################
####State Capacity, Insurgency, and Civil War - Pred. Probs. Country ####
#########################################################################
library(MASS) 
library(pscl)
library(foreign)
library(Hmisc)
library(rgdal)
library(RArcInfo)
library(stargazer)
library(mvtnorm)
library(readstata13)
library(pROC)
library(robust)
library(cvTools)
library(boot)
library(maps)
library(mapdata)
library(doBy)
library(raster)
library(survival)
library(spduration)
library(DAMisc)
library(plyr)
library(DataCombine)
library(Zelig)

########################
###Create The Dataset###
########################
###Importing the main data
#set working directory to location of dta file
setwd("~/Google Drive/NSLC/Final Dofiles")
#Import country leve data
main.data.ct <- read.dta("ntlcount.dta")

##########################
###CW Onset replication###
##########################
#Model 1
log.1.dat <- na.omit(main.data.ct[,c("onset", "lnNL_country_sum", "warl", "lgdpenl1", "lpopl1", "lmtnest", "ncontig", "Oil", "nwstate", "instab", "polity2l", "ethfrac", "relfrac", "second")])
FL.1 <- zelig(onset ~ lnNL_country_sum + warl + lgdpenl1 + lpopl1 + lmtnest + ncontig + Oil + 
              nwstate + instab + polity2l + ethfrac + relfrac,
            data=log.1.dat, model = "logit")
##PPs
range <- seq(min(log.1.dat$lnNL_country_sum), max(log.1.dat$lnNL_country_sum), 0.01)
x.out <- setx(FL.1, lnNL_country_sum = range)
sout <- sim(FL.1, x=x.out)
pdf("pp1.pdf")
ci.plot(sout, ylab="Predicted probability of conflict onset", xlab="Range of nighttime light (country level, in natural log)")
dev.off()

#Model 2
log.1.dat.eth <- log.1.dat[ which(log.1.dat$second > 0.04999),]
FL.2 <- zelig(onset ~ lnNL_country_sum + warl + lgdpenl1 + lpopl1 + lmtnest + ncontig + Oil + 
              nwstate + instab + polity2l + ethfrac + relfrac,
            data=log.1.dat.eth, model = "logit")
summary(FL.2)
##PPs
range <- seq(min(log.1.dat.eth$lnNL_country_sum), max(log.1.dat.eth$lnNL_country_sum), 0.01)
x.out <- setx(FL.2, lnNL_country_sum = range)
sout <- sim(FL.2, x=x.out)
pdf("pp2.pdf")
ci.plot(sout, ylab="Predicted probability of conflict onset", xlab="Range of nighttime light (country level, in natural log)")
dev.off()

#Model 3
log.2.dat <- na.omit(main.data.ct[,c("onset", "lnNL_country_sum", "warl", "lgdpenl1", "lpopl1", "lmtnest", "ncontig", "Oil", "nwstate", "instab", "anocl", "deml", "ethfrac", "relfrac", "second")])

FL.3 <- zelig(onset ~ lnNL_country_sum + warl + lgdpenl1 + lpopl1 + lmtnest + ncontig + Oil + 
              nwstate + instab + ethfrac + relfrac + anocl + deml,
            data=log.2.dat, model = "logit")
summary(FL.3)
##PPs
range <- seq(min(log.2.dat$lnNL_country_sum), max(log.2.dat$lnNL_country_sum), 0.01)
x.out <- setx(FL.3, lnNL_country_sum = range)
sout <- sim(FL.3, x=x.out)
pdf("pp3.pdf")
ci.plot(sout, ylab="Predicted probability of conflict onset", xlab="Range of nighttime light (country level, in natural log)")
dev.off()

